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Parareal convergence for 2D unsteady 
flow around a cylinder 
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Abstract In this technical report we study the convergence of Parareal 
for 2D incompressible flow around a cylinder for different viscosities. Two 
methods are used as fine integrator: backward Euler and a fractional step 
method. It is found that Parareal converges better for the implicit Euler, likely 
because it under-resolves the fine-scale dynamics as a result of numerical 
diffusion. 
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1 Introduction 

The potential of parallel-in-time integration methods to increase the degree 
of concurrency in the numerical solution of time-dependent partial differ¬ 
ential equations has been widely acknowledged, e.g. in the report Applied 
Mathematics Research for Exascale Computing by Dongarra et al. [2014]. A 
variety of different methods exists, see e.g. the review by Gander [2015], and 
principle efficiency of parallel-in-time integration in large- and extreme-scale 
parallel computations has been demonstrated, e.g. in Speck et al. [2012], 
Ruprecht et al. [2013], and Gander and Neumueller [2014]. 

Many problems in computational fluid dynamics require massive compu¬ 
tational capacities and suffer from long solution times. Exploring the poten¬ 
tial of parallel-in-time methods to speed up such simulations can therefore 
be a beneficial endeavour. The performance of the parallel-in-time integra¬ 
tion method Parareal, introduced by Lions et al. [2001], when applied to the 
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Navier-Stokes equations has been a topic of research since shortly after its in¬ 
troduction. First studies have been conducted by Trindade and Pereira [2004, 
2006], including reports of speedup for an MPI implementation; laminar flow 
around a cylinder is used as benchmark problem and it is shown that Parareal 
can correctly reproduce the Nusselt number. Fischer et al. [2005] investigate 
Parareal for spatial discretisations based on finite and spectral element meth¬ 
ods, and discuss using fewer spatial degrees-of-freedom for the coarse integra¬ 
tor. The performance of Parareal for simulations of non-Newtonian fluids has 
been investigated by Celledoni and Kvamsdal [2009]. Finally, parallel scaling 
of Parareal for 3D unsteady flow is investigated by Croce et al. [2014] on up 
to 2048 cores. 

Based on predictions from linear stability analysis by Gander and Vande- 
walle [2007], it has been shown by Steiner et al. [2015]) that convergence of 
Parareal deteriorates as the Reynolds number increases. However, the stud¬ 
ies only analysed a rather simple driven cavity problem, which eventually 
approaches a steady-state and thus may underestimate the problem because 
of weak transient dynamics towards the end of the simulation. In this report, 
we continue this investigation for a different, more complex benchmark in¬ 
volving unsteady flow around a cylinder. It was introduced by Schafer et al. 
[1996] as the case 2D-3 and further analyzed by e.g. John [2004]. Eventu¬ 
ally, the here presented benchmarks will be extended to a comprehensive 
exploration of Parareal’s performance for 3D flow, including a study of the 
influence of spatial resolution. 


2 Parareal and model problem 
2.1 Parareal 

Parareal parallelises the solution of initial value problems 

ut = u{0) = uo, <G[0,T], (1) 

by decomposing the time domain [0,r] into time slices [tj-i,tj], j = 
l,...,lVpr with Npy equal to the number of processing units. It then iter¬ 
ates between two time integration methods: a coarse integrator C used to 
serially propagate corrections, which has to be computationally cheap, and 
an accurate integrator run in parallel. The Parareal iteration reads 

') + - Ciu’l), J = I,... , IVpr, (2) 

with k being the iteration index and Uj « u{tj). Note how the computa¬ 
tionally expensive computation of the fine method can be done concurrently 
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for all time slices. A detailed presentation including a theoretical model for 
projected speedup is given e.g. by Minion [2010]. 


2.2 Model problem 


As model problem, we consider the Navier-Stokes equations 

(dt + UV)U =-—+ lyV'^U (3) 

Q 

for an incompressible fluid, i.e. for a fluid with Vt/ = 0, at density g '■= 
1 (kg/m^). We focus on the benchmark problem defined in Schafer et al. 
[1996] as 2D-3, which is for unsteady flow around a cylinder in two dimensions, 
i.e. in 2D (see also John [2004]). We make use of the definitions 

V:={d,,dyf, U ■.= U{t,x,y), U:=iu,vf. (4) 


The Reynolds number for a cylinder with diameter d (the reference length) 
located inside a square cuboid with longest edge along the x-coordinate is 

- d 

Rd •— rZ-in 5 (5) 

V 

where the reference velocity Uin is chosen to be the mean velocity of inflow 
in a;-direction and n the kinematic viscosity. Notice that Rd can be time 
dependent, as is the case for the problem considered here. 

In the 2D case the mean velocity is 


u{t) := 


u(t,0,y) 
0 h 


dy 





v{t) := 0, 


( 6 ) 


where h := 0.41 (m). In the 2D-3 case we have the inflow velocity 


Min(f, 0,1/) := 4Mf„sin 


y{h - y) 

'■ 


flt,0,y) := 0, 


for which Equation (6) is valid. We choose Mf„ := 3/2 (m/s) so that 


Rd{t) 


/TT \ 

d 

d] 


- e 

0, - 

V8 / 

V 

V 


(7) 

( 8 ) 


Thus, setting d := O.I (m), it follows for n e {0.1,0.01,0.001} (m/s) that 


AS e {1,10,100} 


( 9 ) 


defines the maximum-over-time Reynolds number for the chosen m. 
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In this report, the primary goal is to outline the performance of Parareal 
for the 2D-3 benchmark problem for the three mentioned viscosities, i.e. 
ranges of Reynolds numbers. 


2.3 Implementation details 


The governing equations were implemented using Q2-Q1 finite elements in 
the UG4 software toolbox (see Vogel et al. [2013], Vogel [2014]). For the paral¬ 
lelisation in time via Parareal, we used the library Lib4PrM, which was first 
applied in Kreienbuehl et al. [2015].^ 


3 Results 

Instead of measuring convergence of Parareal by comparing the discretisation 
error with the defect as discussed e.g. by Arteaga et al. [2015], we focus here 
on how well Parareal reproduces important characteristic numbers of the 
dynamics, namely the drag coefficient Cdr, the lift coefficient Cu and the 
pressure difference Ap between the front and end point of the cylinder over 
time. These parameters are respectively defined as follows: 

2Fa 2Fv 

Cdr ■= C'li := •= Pfr - Pen, (10) 

Qufff gufff 

where Fdr is the drag force and F\i the lift force, and pn together with pen 
define the pressure at the front and end of the cylinder. Again, we assume 
here that the density is 1 (kg/m^) and set t € [0, 8] (s) as time domain. 

For n = 0.001 (m/s), Schafer et al. [1996] report on a maximum-over- 
time drag coefficient of « 2.9500 ± 0.0200, a maximum-over-time lift 
coefficient of « 0.4800±0.0100, and a final pressure difference at t = 8 (s) 
of A^ « -0.1100 ± 0.0050 (kg/s2). 


3.1 Numerical setup 

We use Npj. e {2,4,8,16} processors without parallelization in space and 
with 13,212 spatial degrees-of-freedom for both the fine and coarse level. For 
each n, we consider the following two Parareal solvers “S” comprised of a 
coarse Mco and fine Mg serial time integration method as well as number of 


^ It can be obtained by cloning the Git repository https://scm.ti-edu.ch/repogit/ 
lib4prm. 
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time steps Nco and -/Vg: 


(MeciVeo) X (Mfi,iVfi) = (IE, 16) X (IE,32), (Si) 

(Meo,iVeo) X (Mfi, Vfi) = (IE, 16) X (FS,32), (S 2 ) 


where “IE” stands for implicit Euler (first-order) and “FS” for fractional 
step (second-order). Errors in the three physical quantities discussed above 
are measured by 


Eph ■= 




-^phl 


( 11 ) 


where the parameter “ph” is in {dr,li,p} for drag and lift coefficient or pressure 
difference. We use the / 2 -norm over the solutions at the end of all time-slices 


Nr, 


||2 

Ilti 


Nr. 




( 12 ) 


weighted by the time slice length 8/Npr. 


3.2 Problem dynamics 

Figure 1 shows the flow field at f = 5.25 (s) for the three different viscosities. 
As viscosity decreases, the maximum Reynolds number increases and the flow 
becomes more turbulent. While for v = 0.1 (m/s) and v = 0.01 (m/s) the 
flow is essentially laminar, smaller vortices start to form behind the cylinder 
for V = 0.001 (m/s). The resulting evolution over time of the three charac¬ 
teristic numbers for v = 0.1 (m/s) using serial time integration is shown in 
Figure 2. Both fine integrators Ai and S 2 produce essentially identical pro- 
Hles and their profiles closely match the one generated by the corresponding 
reference simulation (not shown). Figure 3 shows the same profiles for the 
simulation with = 0.01 (m/s). Again, both Si and S 2 produce profiles 
that match and agree with the results from the corresponding reference (not 
shown). Lastly, Figure 4 shows three profiles for i' = 0.001 (m/s), one for the 
reference simulation, one for Si and one for S 2 - While the drag coefficient 
and pressure difference agree across all three configurations. Si produces a 
lift coefficient profile that is distinctly different from the corresponding ref¬ 
erence and 5*2. The relatively high numerical diffusion of Si in combination 
with a rather low spatial resolution probably prevents from correctly cap¬ 
turing the more turbulent dynamics in this case. In contrast, although S 2 
fails to fully reproduce the frequency of oscillations, it still achieves a quali¬ 
tatively correct representation of the dynamics of the corresponding reference 
simulation. 
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Fig. 1: Flow field at t = 5.25 (s) for three different viscosities: v = 0.1 (m/s) 
at the top, V = 0.01 (m/s) in the middle, and v = 0.001 (m/s) at the bottom. 


3.3 Convergence of Parareal 

Here, we analyse how accurately Parareal reproduces the three character¬ 
istic values studied above. Figure 5 shows the defect or error according to 
Equation (11) in the characteristic values accumulated over all time slices 
versus the number of iterations. Here, defect refers to the difference between 
the solution computed by Parareal and the solution computed by running 
the fine integrator serially. For Si^ the error for all three quantities, i.e. drag 
coefficient, lift coefficient and pressure difference, quickly goes to zero, that 
is Parareal rapidly produces values identical to ones obtained from the serial 
simulation. As the number of time slices is increased, convergence becomes 
slower but the increase is not drastic: for = 16 after seven iterations all 
three characteristic values have converged up to round-off error. In a produc¬ 
tion run, where the main goal is to push the defect from Parareal below the 
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Fig. 2: Drag and lift coefficient, and pressure difference for ly = 0.1 (m/s) for 
Si (top) and S 2 (bottom). The curves match the values from the reference 
simulation (not shown). 


discretisation error (see the discussion in Arteaga et al. [2015]), significantly 
fewer iterations will likely suffice. 

Decreasing viscosity and thus increasing the Reynolds number range does 
negatively affect convergence. For 1 / = 0.001 (m/s), the simulation with Np^ = 
16 time slices already requires 13 iterations to converge up to round-off error. 
Depending on the desired accuracy, speedup is still possible here but parallel 
efficiency will likely be lower than in the more laminar case. Since fails 
to resolve the full dynamics of the problem, its convergence behaviour is 
probably not representative of the actual physical dynamics. 

This is supported by the fact that for S 2 and v = 0.001 (m/s), Parareal 
essentially no longer converges. Since S 2 does resolve the turbulent dynamics 
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Fig. 3: Drag and lift coefficient, and pressure difference for z/ = 0.01 (m/s) 
for reference (top), Si (middle) and S 2 (bottom). 


at least partially, in contrast to 5i, this suggests that the good convergence 
of Parareal for Si and z/ = 0.001 (m/s) is an artefact produced by excessive 
numerical diffusion. The dynamics of the numerical solution are more lam¬ 
inar than they should be, leading to an unrealistic convergence behaviour 
of Parareeal. Supposedly, when using Si on a significantly finer spatial and 
temporal mesh, a similar deterioration of convergence would be observed, as 
the numerical solution better resolves the turbulent features of the flow. 

Interestingly, this difference between 5*1 and S 2 can already be seen for 
z/ = 0.1 (m/s), where the physical dynamics are still quite laminar as well. 
Although Parareal for S 2 does converge, particularly for larger numbers of 
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Fig. 4: Drag and lift coefficient, and pressure difference for v = 0.001 (m/s) 
for reference (top), Si (middle) and S 2 (bottom). 
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Fig. 5: Fine solver convergence errors for drag and lift coefficient, and pressure 
difference for = 0.1 (m/s) and Si (top) and S 2 (bottom). 


time slices, its rate of convergence is lower than for Si. The benefit of using an 
integrator with damping properties as coarse integrator has been pointed out 
before by Bal [2005] but apparently numerical diffusion from the fine method 
does help Parareal convergence, too. Since here the fine method realistically 
represents the flow features, using a diffusive integrator as fine method in 
Parareal when simulating turbulent flow could be an easy way to obtain 
decent convergence. 
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Fig. 6: Fine solver convergence errors for drag and lift coefficient, and pressure 
difference for = 0.01 (m/s) for (top) and S 2 (bottom). 


4 Conclusions 

This report extends the investigation by Steiner et al. [2015] about how a de¬ 
creasing viscosity in the Navier-Stokes equations affects the convergence be¬ 
haviour of the Parareal parallel-in-time method. An unsteady 2D flow around 
a cylinder is used as a benchmark. Two different configurations are tested, 
one using an implicit Euler as coarse and fine method, the other an implicit 
Euler as coarse but a second-order fractional step integrator as fine method. 
For larger viscosities, both base methods correctly reproduce the evolution of 
characteristic quantities like drag and lift coefficient and pressure difference. 
However, the numerical diffusion from the backward Euler method as fine 
integrator leads to significantly better convergence of Parareal. 
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Fig. 7: Fine solver convergence errors for drag and lift coefficient, and pressure 
difference for = 0.001 (m/s) for Si (top) and S 2 (bottom). 


As viscosity decreases, convergence of Parareal becomes slower similarly 
to the results in Steiner et al. [2015]. However, while for the fractional step 
method Parareal stalls completely, the backward Euler retains reasonable 
convergence even for small viscosities. A comparison of the profiles for char¬ 
acteristic numbers with a reference solution suggests that the good conver¬ 
gence for the backward Euler is likely artificial. The fine integrator alone 
captures the relatively smooth profiles of the drag coefficient and pressure 
difference quite well. However, the high frequency oscillations in the lift co¬ 
efficient, which are clearly seen in serial runs using the fractional step inte¬ 
grator, are not present when using backward Euler. Most likely, the rather 
high numerical diffusion leads to an artificially laminar flow, so that the good 
convergence of Parareal is not representative of the used viscosity parameter. 
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Clearly, when assessing Parareal’s convergence for flow problems, care must 
be taken to ensure that the fine base method correctly resolves the important 
features of the flow. Intrinsic convergence of Parareal alone is not a reliable 
indicator. 

There are several works proposing strategies to stabilise Parareal for 
advection-dominated problems, e.g. by Farhat et al. [2006], Gander and Petcu 
[2008], Ruprecht and Krause [2012] or Chen et al. [2014]. An interesting 
continuation of the work presented here would be to analyse whether these 
strategies improve convergence of Parareal for turbulent flows. 
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